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ABSTRACT 

We model the time variability of ^9,000 spectroscopically eonfirmed quasars in SDSS Stripe 82 
as a damped random walk. Using 2.7 million photometric measurements collected over 10 years, 
we confirm the results of Kelly et al. (2009) and Kozlowski et al. (2010) that this model can 
explain quasar light curves at an impressive fidelity level (0.01-0.02 mag). The damped random 
walk model provides a simple, fast [0(A'') for TV data points], and powerful statistical description 
of quasar light curves by a characteristic time scale (r) and an asymptotic rms variability on 
long time scales (SFoo)- We searched for correlations between these two variability parameters 
and physical parameters such as luminosity and black hole mass, and rest-frame wavelength. 
Our analysis shows SFq^ to increase with decreasing luminosity and rest-frame wavelength as 
observed previously, and without a correlation with redshift. We find a correlation between SFoo 
and black hole mass with a power law index of 0.18±0.03, independent of the anti-correlation 
with luminosity. We find that r increases with increasing wavelength with a power law index of 
0.17, remains nearly constant with redshift and luminosity, and increases with increasing black 
hole mass with power law index of 0.21±0.07. The amplitude of variability is anti-correlated with 
the Eddington ratio, which suggests a scenario where optical fiuctuations are tied to variations in 
the accretion rate. However, we find an additional dependence on luminosity and/or black hole 
mass that cannot be explained by the trend with Eddington ratio. The radio-loudest quasars 
have systematically larger variability amplitudes by about 30%, when corrected for the other 
observed trends, while the distribution of their characteristic time scale is indistinguishable from 
that of the full sample. We do not detect any statistically robust differences in the characteristic 
time scale and variability amplitude between the full sample and the small subsample of quasars 
detected by ROSAT. Our results provide a simple quantitative framework for generating mock 
quasar light curves, such as currently used in LSST image simulations. 
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1. Introduction 

The optical variability of quasars has been rec- 
ognized since they were first identified (Matthews 
& Sandage 1963). Indeed, most quasars are vari- 
able (^90% at the 0.03 mag rms level; Sesar et 
al. 2007), and the variations in brightness are ape- 
riodic and on the order of 20% on time scales 
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of months to years (e.g., Hook et al. 1994; Van- 
den Berk et al. 2004). Furthermore, the smooth 
power spectra suggest a chaotic, or stochastic, ori- 
gin for the variabihty. A range of models have 
been advanced to describe quasar variability, in- 
cluding supernova bursts, microlensing, and ac- 
cretion disk instabilities (Arctxaga ct al. 1997; 
Hawkins 1993; Kawaguchi et al. 1998; Trevese & 
Vagnetti 2002). These models are discussed and 
compared in Hawkins (2007). Reverberation map- 
ping studies (e.g., Peterson et al. 2005) show that 
the broad emission lines respond to continuum 
fluctuations, therefore providing strong evidence 
that the variability is intrinsic to the quasars. A 
number of studies have utilized standard accre- 
tion disk models to demonstrate that the optical- 
UV variability of quasars could be driven by a 
variable accretion rate (e.g., Pereyra et al. 2006; 
Li & Cao 2008; Liu ct al. 2008). Blackburne & 
Kochanek (2010) find evidence in the light curves 
of microlensed quasars that the optical variability 
is caused by a change in the effective area of the 
accretion disk. 

Recently, Kelly et al. (2009, hereafter KBS09) 
proposed a model where the optical variability 
is described by a damped random walk (a self- 
correcting term added to a random walk model 
that acts to push any deviations back towards the 
mean value). They proposed that the variability 
time scale might be identified with the thermal 
time scale of accretion disks, as also proposed by 
Collier & Peterson (2001). A thermal origin of 
the variability would explain why quasars become 
bluer as they brighten (e.g., Givcon et al. 1999; 
Trevese et al. 2001; Geha et al. 2003). 

Although the physical causes have yet to be 
proven, it has been established by KBS09 and 
Kozlowski ct al. (2010a, hereafter KozlO) that a 
damped random walk can statistically explain the 
observed light curves of quasars. Using 100 well- 
sampled single-band light curves compiled from 
the literature, KBS09 show that this stochastic 
process is capable of modeling complex quasar 
light curves at an impressive; fidcility level (0.01- 
0.02 mag). KozlO applied tlie model to the 
OGLE light curves (Udalski et al. 1997; Udal- 
ski et al. 2008) of mid-infrared-selected quasars 
behind the Magellanic Clouds from Kozlowski & 
Kochanek (2009). Their analysis shows that this 
stochastic model is robust enough to efficiently 



select quasars from other variable sources (see 
Schmidt et al. 2010 for a different method of se- 
lecting quasars based on variability). The model 
has only three free parameters: the mean value 
of the light curve, the driving amplitude of the 
stochastic process, and the damping time scale. 
The predictions are only statistical, and the ran- 
dom nature reflects our uncertainty about the 
details of the physical processes. 

Instead of applying a model to observed light 
curves for individual quasars, numerous stud- 
ies have looked at the ensemble variability of 
quasars, particularly in samples where individual 
light curves are not available. Significant progress 
in the description of quasar variability has been 
made by employing the Sloan Digital Sky Survey 
(SDSS) data (Vanden Berk et al. 2004, hereafter 
VB04; Ivezic et al. 2004, hereafter 104; de Vries et 
al. 2005; Wilhite et al. 2005, 2006, 2008; Sesar et 
al. 2006). For example, the size and quality of the 
sample analyzed by VB04 (two-epoch photometry 
for 25.000 spcctroscopically confirmed quasars) al- 
lowed them to constrain how quasar variability in 
the rest frame optical/UV regime depends upon 
rest-frame time lag (up to ~2 years), luminosity, 
rest wavelength, redshift, the presence of radio and 
X-ray emission, and the presence of broad absorp- 
tion line systems. Using repeated SDSS photomet- 
ric observations, Wilhite et al. (2008) confirmed 
the result of Wold et al. (2007) that variability 
is correlated with black hole mass, and show that 
this is independent of the anti-correlation between 
variability and luminosity established by many 
studies. This led them to suggest that the ampli- 
tude of variability may be driven by the quasar's 
Eddington ratio, implying differences in accretion 
rate. 

These studies typically quantify the observed 
optical variability of qiiasars using a structure 
function (SF) analysis (see also Hughes et al. 1992; 
Collier & Peterson 2001; Baii(;r ct al. 2009; 
Kozlowski et al. 2010b), where the SF is the 
root-mean-square (rms) magnitude difference as 
a function of the time lag (At) between mea- 
surements. This autocorrelation-like function is 
less sensitive to aliasing and other time-sampling 
problems than a power spectral distribution. By 
studying the magnitude difference distribution for 
appropriately chosen subsamples with fixed val- 
ues of absolute i-band magnitude {Mi), rest-frame 
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time lag (Atf/p, in days), and wavelength {Xrf, 
in A), the mean dependence of the SF on these 
quantities was inferred by 104 to be 

/At \ *^ 

SF^,aei = A[l + BM,]i—^] mag, (1) 

\ ^RF / 

with A ^ 1.00 ± 0.03, B = 0.024 ± 0.04, and 
C = 0.30 ±0.05. A qualitatively similar result was 
obtained by VB04. Kozlowski et al. (2010b), in 
the first large study of the mid-IR structure func- 
tions of quasars, also found lower variability for 
higher luminosities and longer wavelengths, but 
the temporal slope of the ensemble structure func- 
tions were significantly steeper than in the optical. 
In addition, there is evidence for a turnover in the 
SF on long time lags (104; Rengstorf et al. 2006; 
Wold et al. 2007). Studies by de Vries et al. (2005) 
and Sesar et al. (2006) using SDSS combined with 
earlier Palomar Observatory Sky Survey measure- 
ments for 40,000 SDSS quasars constrained quasar 
continuum variability on time scales of 10 to 50 
years in the observer's frame. They report that 
the characteristic time scale, which in this context 
is the time lag above which the SF flattens to a 
constant value, is of order 1 year in the quasar 
rest frame. Using a shot-noise light curve model, 
de Vries et al. (2005) found evidence for multi- 
ple variability time scales in long-term ensemble 
variability measurements, while Collier & Peter- 
son (2001) found a wide range of different time 
scales in their analysis of individual light curves, 
and even evidence for multiple time scales in a 
single active galactic nucleus (AGN). 

These analyses of ensemble variability are based 
on a fundamental assumption that photometric 
observations at two epochs for a large number of 
quasars will reveal the same statistical properties 
as well-sampled light curves for individual objects. 
This assumption has been tested by MacLeod et 
al. (2008) using light curves for spectroscopically 
confirmed quasars observed roughly 50 times over 
8 years in SDSS Stripe 82 (S82). They found 
that while the mean SF for individual sources 
is consistent with Eq. [TJ the contribution of the 
mean trends to the observed dispersion in variabil- 
ity properties is minor compared to an intrinsic 
stochasticity of unknown origin. Further inves- 
tigation of this stochastic behavior is one of the 
main goals of this study. 

In order to better understand the relationship 



between the two types of data analyses (individual 
versus ensemble quasar variability), and to begin 
linking to physical models, we apply the damped 
random walk model to the ugriz light curves 
of ~9,000 spectroscopically confirmed SDSS S82 
quasars. This large sample greatly benefits from 
the robust, accurate, five-band SDSS photometry. 
We estimate the variability parameters following 
KozlO, who demonstrated that their approach is 
more statistically powerful than the forecasting 
methods used by KBS09. We also note that both 
the KozlO and KBS09 approaches are much faster 
than that used by Schmidt et al. (2010), requir- 
ing only 0{N) rather than 0{N'^) operations to 
determine the model parameters for a light curve 
with N data points. 

In Section [21 we describe the model, define our 
variability parameters, and demonstrate their re- 
lationship to those utilized in previous studies. In 
Section [21 we introduce the S82 data set and out- 
line our initial light curve selection. In Section [4l 
we present the best-fit variability parameters for 
our final sample of light curves and estimate their 
scatter due to the limited time sampling of SDSS. 
We also estimate the sensitivity of our results to 
variations in the slope of the model power spec- 
trum on long time scales. In Section[5l we describe 
the relationship between the long-term variabil- 
ity parameters and physical parameters such as 
wavelength, absolute magnitude, black hole mass, 
and Eddington ratio. Using these results, we also 
provide a prescription for simulating mock quasar 
light curves. In Section [6l we explore the variabil- 
ity properties of subsamples detected at radio and 
X-ray wavelengths. Finally, we summarize our re- 
sults in Section [3 

2. Methodology 

We model the time variability of quasars as a 
stochastic process described by the exponential co- 
variance matrix 

^a^eM-\U-tj\/T) (2) 

between times ti and tj. As detailed by KBS09 
and KozlO, this corresponds to a damped random 
walk (more specifically, an Ornstein-Uhlenbeck 
process) with a damping time scale t, also called 
the characteristic time scale, and a long-term stan- 
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dard deviation of variabilitj0 a. Following KozlO, 
we model the light curves and estimate the param- 
eters and their uncertainties using the method of 
Press et al. (1992), its generalization in Rybicki & 
Press (1992), and the fast computational imple- 
mentation described in Rybicki & Press (1994). 
KozlO demonstrate that this approach is more 
statistically powerful than the forecasting meth- 
ods used by KBS09, while still having computa- 
tion times scaling linearly with the number of data 
points. 

2.1. Structure Function for the Damped 
Random Walk Model 

For our analysis, we express the long-term vari- 
ability in terms of the structure function (SF) 
in order to relate to our previous studies and 
to those of two-epoch samples. The first or- 
der SF is SF^{At) = 2cr2[i _ ACF{At)] (e.g., 
Hughes et al. 1992), where the autocorrelation 
function for a damped random walk is ACF{At) = 
exp(-|At|/r) (Eq. A6 in KBS09). This ACF re- 
sults in the structure functiorH: 

5^^(At) = 5^00(1-6-1^^*1/^)1/2. (3) 

Asymptotic values of the SF at large and small At 
are 

SF{At »t) = SFoc = V2ct 
SFiAt«r)=a^^SF^^ (4) 

The form SF cx \At\^ is equivalent to a power 
spectral distribution PSD cx where a = 

—2/3—1 (see Appendix of Bauer et al. 2009; 
KBS09). The SF at smaU time lags is there- 
fore equivalent to a power spectral distribution 
PSD cx /~2. We adopt SFoo and r as our two 
main variability model parameters. 

2.2. Model Light Curves 

Equipped with a statistical description of 
quasar variability, we generate well-sampled light 



^The (T used here is related to the a used in KBS09 (urbs) 
and the parameter a used in KozlO as crj^^gg = a = 
cry2/V. 

2The functional form of SF(Ai) fit to the long-term SDSS- 
POSS data in Sesar et al. (2006), see their Eq. 5, is similar 
but not identical to the functional form given by Eq. [3l 



curves in order to 1) demonstrate the relation- 
ship between our variability parameters and the 
traditional SF analyses of many previous works, 
and 2) to estimate the systematic effects that the 
sampling rate and light curve length have on the 
fitted parameters. The latter is especially impor- 
tant because the S82 data are fairly sparse. As 
shown below and in Section |4?2| these indeed have 
a large impact. 

A light curve is generated using only three input 
parameters: SFoo, t, and the mean value of the 
light curve, /i. The magnitude X(t) at a given 
timestep At from a previous value X(t — At) is 
drawn from a normal distribution with a mean 
and variance given by 

E{X{t)\X{t - At)) = e-'^'/"X(f - At) + 

/I —At/T\ 

m(1 - e ' ) 
Var{X{t)\X{t - At)) = 0.5{SF^)^{1 - e""^'/") (5) 

(Eqs. A4 and A5 in KBS09). The asymptotic 
variance of the time series is then 0.5(5^00)^. 
The top panel of Figure [T] shows a segment of a 
well-sampled light curve generated using r = 20 
days, SFoo = 0.14 mag, and a time sampling of 
0.1(At/r). The structure function, SF{At), is 
computed by collecting the differences in magni- 
tude for all points in the light curve separated by a 
given time lag. At. The distribution of magnitude 
differences (Am) is Gaussian by construction, and 
the rms of the Am distribution is the SF value for 
that time lag. 

When fitting values for SFoo and r for a given 
light curve, the length of time that it spans plays 
an important role. For example, in Figure [TJ 
the SF computed for the full light curve length 
of 1500r (82 years) is much smoother than that 
computed for three equal sections, each spanning 
27 years. Therefore, the determination of vari- 
ability parameters for S82 quasars will be affected 
by their light curve lengths, which are typically 
lOr. When the light curve is too short, it is 
easy to overestimate r because there is no in- 
formation on the time scale of the break (where 
SF{At) flattens to SFoo)- The model will repro- 
duce the observed variance in the light curve by 
overestimating SFoo — {S Foo)true{T / At)^^'^ , and 
therefore a — SFoo / \/t is the more robustly es- 
timated model parameter when t cannot be well- 
determined. 
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2.3. Comparison with Published Work 

Before interpreting the form of SF{At), we 
summarize the major differences between our S82 
analysis and previous studies based on ensemble 
structure functions. The ensemble structure func- 
tion is computed using only a few observations 
of many quasars, combining all magnitude dif- 
ferences to find SF{At). Using an ensemble ap- 
proach is beneficial because it enables one to con- 
strain the average variability properties when it 
is difficult to constrain such quantities for indi- 
vidual quasars. Indeed, even with well-sampled 
light curves, spurious breaks in the individual SFs 
are common (Emmanoulopoulos, McHardy, & Ut- 
tley 2010). However, in previous works (e.g., 104; 
de Vries et al. 2005), the characteristic time scale 
is defined as the time lag at which the ensemble 
SF{At) flattens to a constant value, and thus it 
represents some complex average over the intrinsic 
r distribution. In contrast, by applying a stochas- 
tic model to the individual S82 light curves, we 
are relatively insensitive to time sampling issues 
(for details, see KBS09), and we obtain a model 
fit for every quasar in each filter described by the 
parameters SFoo and r. 

A power law fit to SF{At) has previously been 
a common way to describe SF{At), and even to 
reject certain classes of models (e.g., Kawaguchi 
et al. 1998; VB04; 104). However, Figure □ shows 
that the best-fit power law index is extremely sen- 
sitive to the fitted range of At/r. For example, for 
0.1 < At/r < 1, SF{At) is well fit by a power law 
with an index of ~0.5, while for 0.15 < At/r < 3, 
we obtain a strongly biased power law index of 0.3. 
Furthermore, each quasar has its own values of r 
and iSFoo, and the ensemble structure function is 
a convolution of the individual structure functions 
with their distribution in parameters, 

dTdSFoo—r^SF{At\r,SF^), (6) 

where SF{At\T, SFoo) (Eq. [S]) is the structure 
function at time At for a quasar with variabil- 
ity parameters t and SFoo- Hence, the ensemble 
structure function is only indirectly related to the 
structure function for any particular quasar, and 
results based on fitting a power law to observed 
ensemble structure functions should be interpreted 
with caution. 



3. The SDSS Stripe 82 Quasar Data Set 

The Sloan Digital Sky Survey (SDSS, York 
et al. 2000) provides homogeneous and deep 
(r < 22.5) photometry in five passbands {ugriz, 
Fukugita et al. 1996; Gunn et al. 1998; Smith et 
al. 2002) accurate to 0.02 mag, of almost 12,000 
deg^ in the Northern galactic cap, and a smaller, 
but deeper, survey of 290 deg^ in the Southern 
galactic hemisphere. For this 290 deg^ area known 
as Stripe 82 (S82), there are on average more 
than 60 available epochs of observations. These 
data were obtained in yearly "seasons" about 2-3 
months long over the last decade, and the cadence 
effectively samples time scales from days to years. 
The light curve lengths are effectively shorter than 
the actual period of the survey because the better- 
sampled supernova observations begin about 5 
years into the survey. Because some observations 
were obtained in non-photometric conditions, im- 
proved calibration techniques have been applied 
to SDSS S82 data by Ivezic et al. (2007) and Sesar 
et al. (2007), and we use their results. For these 
data, photometric zero-point errors are 0.01-0.02 
mag. 

We have compiled a sample of 9,275 spec- 
troscopically confirmed quasars in S82 with re- 
calibrated ugriz light curves (see also Bhatti et 
al. 2010). Most (8,974) of these are in the SDSS 
Data Release 5 (DR5) Quasar Catalog (Schnei- 
der et al. 2007), and the remaining are newly 
confirmed DR7 (Abazajian et al. 2009) quasars. 
Summed over all bands and epochs, the data set 
includes 2.7 million photometric measurements. 
For 41% of the sample, the random photomet- 
ric errors are smaller than 0.03 mag. Only 1% 
have errors > 0.1 mag in g, r, and i, and 2.4% 
have errors exceeding 0.25 mag in u and z fil- 
ters. In MacLeod et al. (2010, in prep.), these light 
curves, as well as a much larger sample of quasars 
with two SDSS epochs selected from 12,000 deg^ 
of the sky, will be made publicly available. We 
adopt the K-corrected i-band absolute magnitudes 
from Schneider et al. (2007), and virial black hole 
masses and bolometric luminosities where avail- 
able from Shen et al. (2008). The Shen et al. 
masses were estimated from emission line widths 
(H/3 for z < 0.7, Mgll for 0.7 < z < 1.9, and CIV 
for z > 1.9). However, we note that at low spec- 
troscopic signal-to-noise, black hole masses tend 
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to be overestimated (Denney et al. 2009). 

3.1. Initial Light Curve Selection 

The damped random walk model was fit to all 
available ugriz light curves for 9, 275 S82 quasars. 
Summed over 5 bands, there are 46,375 best-fit 
values of the characteristic (damping) time scale 
r and long-term structure function SFoc- For fur- 
ther analysis, we select light curves that satisfy the 
following criteria: 

1. First, we remove light curves with fewer than 
10 observations. The top-left panel of Fig- 
ure [2] shows the distribution of the number 
of observations (Nobs) per light curve be- 
fore this cut, which reduces our sample to 
45,814. At a given Nobs, the distribution of 
the ratio of light curve length to r is similar 
to that at all other values of Nobs ■ Therefore, 
any systematic effects of Nobs on derived pa- 
rameter distributions should be small. 

2. We then require that the stochastic model 
must provide a better fit than uncorre- 
lated noise. Following KozlO, we select 
light curves with a likelihood improvement 
over simply broadening the measurement 

errors of ALnoise = In (ibest/inoise) > 2, 

where Lbest is the likelihood of the stochas- 
tic model and Lnoise is that for the noise 
solution where r — ^ 0. The top-right panel 
of Figure [2] shows the distribution of l^Lnoise 
before this cut, which removes 14% of u and 
z light curves, whose photometric errors are 
larger. About 7% of our light curves (82% 
of which are u or z band) are removed in 
this stejl^, reducing our total to 42, 623. 

3. Finally, we remove cases where r is merely 
a lower limit due to the length of the light 
curve. We define ALoo = In (ibest/^oo), 
where Loo is the likelihood that r — >■ oo, in- 
dicating that the light curve length is too 
short to accurately measure r. The bottom- 
left panel shows a peak at ALqo = 0, and 
we exclude these objects by requiring that 
ALoo > 0.05. Most (95%) of the rejected 



A detailed analysis of the use of this variability model to se- 
lect candidate quasars will be presented elsewhere (Brooks 
et al., in prep.). 



light curves have lengths < r; 76% have 
T > 10** days, and 64% have SFoo > 1 mag. 
The latter is due to the fact that as the r 
value becomes long and uncertain, the model 
will necessarily overestimate 5'L'oo in order 
to keep the overall light curve rms fixed (see 
Section [221). 

The rejected light curves tend to be higher 
redshift quasars because stronger time dila- 
tion leads to shorter rest-frame light curve 
lengths, making it increasingly difficult to 
constrain long rest-frame r. This criterion 
removes 22% of our sample, leaving a total 
of 33, 112 values of r and SF^o- Because we 
are limited by the duration of the S82 survey, 
this is a significant loss, and therefore our fi- 
nal T distribution is biased low, but the bias 
should not be significant considering our re- 
sults in Section 221 

The resulting light curves are well-fit by the 
stochastic model, as can be seen from the distribu- 
tion of x^/iVdo/ shown in the bottom-right panel 
in Figure [21 where N^of is the number of degrees 
of freedom. The expected Gaussian distribution 
with rms = 2 /N^of is also shown in the panel, 
where we have averaged over the Ndof distribu- 
tion of the light curves. The observed distribution 
is centered at /Ndof = 1-1. This difference is 
some combination of errors in the estimated er- 
rors, outliers in the light curves, and any poorly 
modeled physics. KozlO noted a similar difference 
in their analysis of OGLE light curves. Only 5% 
of the light curves have /Ndof > 1.5, confirming 
that most quasars are variable (at the AL„oise > 2 
level), and that a damped random walk is a good 
description of quasar variability. 

4. Variability Parameters for Stripe 82 
Quasars 

4.1. Observed Distributions 

We assume the variability to be intrinsic to the 
quasars and convert the time scales to the rest 
frame (dividing by (1 -I- z)) before further analy- 
sis. The long-term structure function should be 
independent of redshift other than through evolu- 
tion in physical parameters, and this view is con- 
firmed by KBS09. Figure [3l shows the distribu- 
tions in SFoo and rest-frame r found for the S82 
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quasars. If we consider only the brighter {i < 19) 
quasars, based on the best-lit mean magnitude, 
the distributions are generally similar but biased 
towards lower asymptotic amplitudes (peaked at 
0.12 mag). The distributions are consistent with 
what was found in KBS09 and KozlO; however, 
in these studies, there are not as many objects 
with runaway (t — > cxd) time scales because of 
the improved time sampling of the OGLE survey 
and many of the light curves in the KBS09 sam- 
ple. Also, the observed-frame r distribution lacks 
many of the short time scales observed in KozlO; 
this is likely due to either the better time sampling 
of the OGLE light curves or unrecognized stellar 
contamination in their sample. 

The bottom panel of Figure |3] shows that the 
best-fit variability parameters are highly corre- 
lated with each other, indicating that quasars with 
larger asymptotic amplitudes of variability also 
have longer characteristic time scales. We fit a 
power law slope of 1.3±0.01 dex/dex for all 33,112 
data points, and this trend persists within each 
ugriz band as well. Note that a correlation be- 
tween T and SFoc is expected even if r is indepen- 
dent of the driving amplitude of short-term varia- 
tions, a = SFoo/\/t. Since the power law slope of 
1.3 is steeper than that expected if r and a are in- 
dependent (^0.2, see below), the time scale must 
also be correlated with the amplitude of short- 
term variability. We confirm that this correlation 
is intrinsic, rather than an artifact of short light 
curve lengths, in Section |4?2] 

The approximate values for r and SFoo can of- 
ten be guessed from light curves by simple visual 
inspection. Figure |4] shows representative light 
curves from several regions otr— SFoo space, indi- 
cated by stars in Figure [3] The weighted average 
of all model light curves consistent with the data is 
also shown, and the "error snake" is the ±la range 
of those light curves about this mean. The top 
panel shows a light curve with a relatively short r - 
this is due to the large amount of variability within 
each season (i.e. each grouping of points), as com- 
pared to that in the second panel, which shows a 
curve with a much larger r. The third panel shows 
a light curve with a relatively low SFoo , while the 
bottom panel shows one with a larger SFoo, as 
can be seen by the larger difference in median r 
band magnitude between about 10 and 11 years. 
Note the outlier in the third panel: outliers such as 



these are generally responsible for higher values of 
X^/^dof- The (data— model) brightness difference 
provides a convenient way to identify outliers, and 
they may be an indicator of variability behavior 
not captured by the damped random walk model. 
However, they may also be caused by occasional 
non-Gaussian photometric errors and thus their 
analysis requires a detailed and careful study (e.g., 
by utilizing control samples of appropriately cho- 
sen nearby non- variable stars). Since the model 
provides satisfactory fits for the overwhelming ma- 
jority of objects (see Section [5TT|) . we do not fur- 
ther investigate such outliers in this work. We 
have also searched for periodic signals in observed 
light curves and did not find any convincing cases. 
This analysis is summarized in the Appendix. 

4.2. The Effect of S82 Time Sampling and 
Estimate of Fitting Errors 

There are three contributions to the scatter in 
the best-fit variability parameters t and SFoo- 
First, the fitting errors, including those due to 
insufficient time sampling and light curve length, 
will introduce some scatter. Second, trends with 
physical parameters (see Section [5]) result in a cer- 
tain distribution width. There may also be scat- 
ter due to other sources of variability that are not 
captured by the model, such as flares, or other ac- 
tivity related to radio emission, for example. Here, 
we carry out two tests in order to understand how 
the best-fit parameters are affected by the limited 
data sampling of S82. 

Since the correlation between r and SFoo seen 
in Figure [3] might be expected if the light curves 
are not sampled over sufficiently long periods of 
time (see Section \2.2^ . we first test whether the 
correlation is real or simply an artifact. We gen- 
erated light curves as described in Section [22] us- 
ing the time sampling and photometric uncertain- 
ties of the S82 data. For each object, the in- 
put parameters Tin and SF{oo)in are randomly 
drawn from a uniform distribution in log r and 
log 5*^00 limited by the dotted rectangle in the 
left panel of Figure [5] These artificial light curves 
are then fit to obtain Tout and SF{oo)out, and 
the resulting distribution is shown by the con- 
tours in Figure [5] The open circles show data 
points that do not satisfy ALnoise > 2 (3% of 
all input points) - these are concentrated at small 
Tout and SF{oo)out- The closed circles show those 
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that do not satisfy ALqo > 0.05 (21% of all input 
points) - these have been smeared to large val- 
ues of Tout and SF{oo)out- After omitting points 
with ALnoise < 2 and ALqo < 0.05, the distribu- 
tion of the output estimates is similar to the input 
distribution and shows no strong correlation be- 
tween Tout and SF{oo)out- This suggests that the 
correlation in Figure [3] is largely real, and not an 
artifact of sampling and fitting procedures. These 
results also justify the selection cuts outlined in 
Section 13.11 The right panel of Figure [5] shows 
the expected correlation that the larger the over- 
estimate of T, the larger the overestimate of SFoo , 
with a slope between them following that expected 
by Eq. m (for At << r). We repeated the test 
using a uniform input distribution in log a (the 
driving amplitude of short-term variations) rather 
than in logS'Foo, where a — SFoo/\/t. In this 
case, we find the output t and SFoo are corre- 
lated with a power law slope of 0.18 ± 0.01. Since 
this slope is smaller than that for the observed 
distribution in Figure [3l the a and t must be in- 
trinsically correlated for the S82 sample as well. 

For our second test, we used the best-fit r, 
SFoo, and fi for the S82 sample to generate new 
light curves with the same time sampling and 
photometric uncertainties as the S82 data. By 
comparing the output and input parameter dis- 
tributions, we can estimate how much the intrin- 
sic stochasticity and the time sampling issues af- 
fect the results. Figure [5] compares an observed 
and "regenerated" light curve, where the differ- 
ences are due to the stochastic nature of the pro- 
cess. The fit parameters for the "regenerated" 
light curve can be very different because of how the 
particular realization is affected by the time sam- 
pling, as illustrated in Figure [6] Figure [7] shows 
the ratio of output to input distributions for both 
r and SFqo. The input distributions normalized 
by their median values are also shown to illustrate 
their dynamic range. These two distributions are 
compared to each other in order to estimate the 
effect of fitting errors (the ratio of output to in- 
put should be a delta function centered at 1 for 
perfect time sampling). The bottom-right panel 
shows that the correlation between r and SFoo 
becomes slightly weaker than that seen in Fig- 
ure [3] Based on these results, we conclude that 
the uncertainties due to sparse sampling and lim- 
ited lengths of the S82 light curves can account 



for 71% of the spread in SFoo and 57% of that in 
T. As shown in Section [2?2| very long light curve 
lengths are needed to estimate accurate time scales 
and asymptotic amplitudes. Nevertheless, the ob- 
served distributions indicate that the underlying 
intrinsic distributions of r and SFoo have finite 
widths that are similar to the observed widths. 

4.3. Relationship between the Individual 
SFs and the Ensemble SF 

The distribution of r is important to consider 
when relating the ensemble SF, such as those de- 
termined using two-epoch datasets, to the SFs for 
individual light curves. The analysis based on two- 
epoch data measures the mean value of the SF, 
but provides no information about the SF vari- 
ance among individual objects. To measure the 
latter, individual light curves must be available. 

MacLeod et al. (2008) analyzed individual light 
curves for S82 quasars in order to test the com- 
mon assumption that photometric observations at 
two epochs for a large number of quasars will re- 
veal the same statistical properties seen in light 
curves for individual objects. They found that 
the dependence of the mean SF computed using 
SFs for individual light curves on luminosity, rest- 
frame wavelength and time lag is indeed qualita- 
tively and quantitatively similar to that derived 
from two-epoch observations of a much larger sam- 
ple. However, they also found that the scatter in 
the light-curve based SFs for fixed values of Mi, 
Xrf and Atup is very large, and in fact, similar 
to the scatter for the whole sample (see Figure [8|). 
This large scatter was attributed to an intrinsic 
stochasticity of unknown origin. The new model- 
based analysis discussed here allows us to explain 
this puzzling result as a consequence of the finite 
width of the r distribution. 

An important piece of information missing from 
the MacLeod et al. (2008) analysis is the existence 
of a characteristic time scale r. In MacLeod et 
al. (2008), the individual structure functions (SF) 
were computed following the standard approach 
with a fixed observer's frame time lag of 1 year 
(and with rest-frame time lags spanning 100 to 
300 days). However, Eq. 2] indicates that the SF 
at small time lags should be proportional to l/y^, 
indicating that the observed structure functions 
will vary between quasars even if they have similar 
SFoo- Indeed, the distribution of ■y/niedian(r)/T 
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has an almost identical shape and width as that 
of SF/SFmodei (see Figure H]). We therefore con- 
clude that variations in r are responsible for most 
of the scatter in SF{At << r) for quasars with 
similar luminosity, rest wavelength, and time lag. 
In addition, it is likely that the r distribution is 
responsible for the exponential tails of the mag- 
nitude difference distribution for quasars reported 
by 104 and Sesar et al. (2006; see MacLeod et al. 
2010, in prep., for further discussion). Therefore, 
the published structure function results based on 
two-epoch datasets can only be interpreted in the 
context of Eq. [6l 

4.4. Dependence on the Underlying PSD 

The analysis throughout this paper assumes 
that the variability is described by a damped ran- 
dom walk, which has a power spectral distribution 
(PSD) described by PSD cx at frequencies 
/ > (27rr)~^, flattening to a constant at lower 
frequencies. In this Section, we investigate the 
sensitivity of our resulting parameter distributions 
to the possibility that the low frequency part of 
the PSD is not flat. For example, the X-ray vari- 
ability of Seyfert galaxies is well-described by a 
broken power law with a slope of —2 at high fre- 
quencies, breaking to a shallower slope (—1) at low 
frequencies (Arevalo et al. 2008a). Since the opti- 
cal and X-ray variations are correlated (Arevalo et 
al. 2008b), it is plausible that the optical variabil- 
ity might have a similar underlying PSD. On the 
other hand, as reviewed by McHardy (2010), the 
optical and X-ray fluctuations are only correlated 
on time scales which are shorter than the typi- 
cal characteristic optical time scale (i.e., on time 
scales corresponding to the 1/P part of the opti- 
cal power spectra) . If the damped random walk is 
a good description for the optical variability, then 
we might expect that the optical and X-ray fluc- 
tuations are no longer correlated on time scales 
longer than r, as the optical fluctuations resemble 
white noise on these time scales. 

We consider three cases for a broken, or bend- 
ing, PSD described by PSD oc /^^ at high fre- 
quencies and PSD oc /" at low frequencies. The 
first case is where a = 0, which is a damped ran- 
dom walk. The second case is where a = — 1, and 
the third is where a = —1.9, which is nearly a 
constant power-law slope. In each case, ~7000 
light curves are simulated using the algorithm 



from Timmer & Koenig (1995) with the chosen 
PSD. For each realization, the break frequency is 
set to {2ttt)~^ and the total rms is fixed, using 
the observed r-band r and rms values for the S82 
quasars that satisfy the selection criteria in Sec- 
tion [3?ll Therefore, the observed distributions in 
Figure [3] (for the r-band) are the "input" values. 
Figure [9] shows 3 example light curves simulated 
using a — ~1 (red), a = (blue), or a = —1.9 
(green). The lines in the last panel show each PSD 
from which the light curves are generated, with the 
black dotted line indicating the break frequency. 
As a check, the PSD was computed for each simu- 
lated light curve, and as seen by the colored dots, 
the shapes match the input. 

Each light curve was simulated over 100 years 
and then truncated to the 30-40 year segment in 
order to account for the additional variability and 
bias in the inferred SFo^ , which may result from a 
red noise leak or an incorrect specification of the 
PSD model (see Uttley et al. 2002). The simu- 
lated light curves were then modeled as a damped 
random walk to obtain the ("output") parameter 
distributions shown in Figures fTOl and fTTl For Fig- 
ure (TUl the simulated observations are spaced ev- 
ery 5 days over 10 years with typical errors of 0.01 
mag, and for Figure [TTl the S82 window function 
is imposed (i.e., all light curves have the S82 time 
sampling and photometric accuracy). The filled 
histograms show the input distributions and are 
the same for both Figures. 

It is clear from comparing the input and out- 
put parameter distributions that r and SFoo map 
simply onto the values of the power-spectral break 
time scale and amplitude. Indeed, the correlation 
between r and SFoo is preserved in the output 
distributions, as seen in the bottom-right panels 
of Figures \W\ and [TT] However, as seen in the 
top- and bottom-left panels, as the PSD slope at 
low / steepens to a = —1.9, the number of "run- 
away" time scales (where t saturates at 10^ days) 
increase due to that fact that r can no longer be 
constrained. Whereas for the S82 data, ~20% of 
light curves were rejected as run-away cases (see 
step 3 of Section [3TT|l . for an a = —1.9 PSD, the 
run- away fraction is 30%. This significant increase 
rules out the a = —1.9 PSD as the correct model, 
as we do not see this large run- away fraction in 
the data. The fractions are similar for the a = — 1 
and a — cases, suggesting that the correct model 
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has — 1 < a < 0. The increasing fraction of run- 
away T with decreasing a ahows one to distinguish 
between each PSD in the well-sampled case, as 
seen in the top-right panel of Figure [lUl In this 
panel, the difference between the Xpdf the best- 
fit damped random walk and that for a r — ;> cxd 
solution, Axlo = X^l^dof ~ xL/Ndof, is shown 
for each input PSD. The total distribution shows 
two peaks; that on the left corresponds to cases 
where the model is able to constrain r. For well- 
sampled light curves (Figure [10]), the case with 
a = (dashed line) yields an overall lower 
and therefore a more likely damped random walk 
fit, as expected. However, for the S82 sampling, it 
is difhcult to distinguish between the red and blue 
lines {a — and —1), and the (input) observed 
distribution is similar to both. Therefore, within 
the limited S82 sampling, we are unable to dis- 
tinguish reliably between a damped random walk 
and a 1// PSD on long time scales. 

5. Dependence of Variability Parameters 
on Luminosity, Wavelength, Redshift, 
and Black Hole Mass 

Next, we discuss correlations between each of 
the variability parameters, t and SFoo, and the 
four available physical parameters: rest-frame 
wavelength (Xrf), redshift (z), absolute i-band 
magnitude (Mi), and black hole mass [Mbh]- 
It is important to fit a multiple regression to 
each of the physical parameters, functions of the 
form T{XiiF , Mi, Mbh , z), because of the correla- 
tions between physical parameters. For example, 
when searching for a correlation between SFoo 
and Mbh, one must take into account that a 
more luminous quasar hosts a more massive black 
hole (e.g., KoUmeier et al. 2006), or else a trend 
with luminosity may be mistaken for a trend with 
black hole mass. A similar example is the well- 
known luminosity-redshift {L~z) degeneracy seen 
in flux-limited samples of quasars. Magnitude 
limits result in the illusion that only the most lu- 
minous quasars are seen at high z, and therefore 
the observed L increases with z independent of 
reality. Having a large sample size helps to alle- 
viate these degeneracies. For example, one can 
look for trends using 2-dimensional grids in any 
two physical parameters of interest. The numbers 
of data points per two-dimensional bin of Mi and 
redshift, and similarly for Mi and Mbh, are shown 



in the bottom panels of Figure [121 The two left 
panels of Figure [T^ show the selection effect that 
quasars at higher redshift must have higher lumi- 
nosities to be included in the survey. From the 
two-dimensional distribution in the top-left panel, 
we can see that this L-z degeneracy is essentially 
independent of Mbh- Furthermore, shorter rest- 
frame wavelengths are probed at higher redshifts 
for two reasons: first, quasars emitting at shorter 
rest wavelengths must be at higher redshifts in 
order to be observed within the ugriz filters, and 
second, quasars at high redshifts are closer to their 
Eddington limit as a result of the L-z degener- 
acy, and possibly cosmological downsizing (e.g., 
Kollmeier et al. 2006). Therefore, any dependence 
of variability on wavelength must be accounted for 
when considering a dependence with redshift (or 
luminosity) . 

When fitting power laws to a large number of 
data points throughout this paper, we fit to the 
median values in each bin of the independent vari- 
able, where all bins have the same number of data 
points N . The errors in the medians are computed 
as 0.93(1QR)/VA^- 1, where IQR is the 25% - 
75% interquartile range, and these constrain our 
formal errors in the power law slopes. 

5.1. Trends with Rest-frame Wavelength 

We start by examining the wavelength depen- 
dence of the variability parameters. Since there 
are multiple bands for each quasar, this depen- 
dence can be determined for individual sources 
for which z, Mi, and Mbh are fixed. The rest- 
frame wavelength is found by dividing the ob- 
served bandpasses (3520, 4800, 6250, 7690, and 
9110 A for ugriz, respectively) by (1 -I- z). We 
fit a power law / oc {\rf /^^^^A)^ to the esti- 
mates of T and SFoo for every quasar observed in 
at least two filters (^^8,000 quasars). The median 
values are B = 0.17±0.02 and -0.479±0.005 for r 
and SFoo , respectively. We searched for significant 
correlations between B and other physical param- 
eters and did not find any. We use these median 
values to correct for the wavelength dependence of 
T and SFoo before searching for their correlations 
with other physical parameters in Sections l5.2[|5.3l 
and [6] 

This method of fixing B before investigating 
other correlations naturally eliminates any degen- 
eracies between rest wavelength and the other 
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physical parameters. This is especially impor- 
tant in the case of the variability amplitude, SFoo- 
Figure [13] shows how the two variability parame- 
ters vary with X^p in each of the ugriz filters. 
The median power law from the above analysis, 
shown as a straight line stretching across all wave- 
lengths, accurately traces the overall trend. How- 
ever, when the data are fit to ensembles of quasars 
in each band separately, the slopes for SFrx, , shown 
by the short lines for each filter, are very different 
(~0.4). This difference is a consequence of the 
correlation between L and z: for a given band, 
shorter Xrf corresponds to higher z, but at higher 
z, quasars have higher L and thus smaller SFoo 
(see below), creating a bias in the inferred slope 
of the wavelength dependence. 

5.2. Trends with Luminosity, Redshift, 
and Black Hole Mass 

In the upper panels of Figure I14[ the median 
values of SF^o and r are shown as a function 
of absolute magnitude {Mi) and redshift. The 
structure function parameters are normalized to 
a fixed rest wavelength of 4000A using the fitted 
power law dependencies of (Ai?F/4000A)^ with 
B = —0.479 and 0.17 for SFqo and t, respectively, 
from Section FS.ll before finding the median in each 
pixel. For SF^o, the anti-correlation with lumi- 
nosity clearly dominates any trend with redshift. 
The T distribution also shows negligible correla- 
tion with redshift. The bottom panels show the 
dependence on Mbh- Using a grid of Mi vs. Mbh 
allows one to search for trends in the variability 
parameters while accounting for the selection ef- 
fect that more luminous galaxies host more mas- 
sive black holes. A positive correlation between 
SFoo and Mbh is apparent, independent of the 
correlation with Mi . This is in agreement with the 
result from Wilhite et al. (2008), who used ensem- 
ble structure functions. The t parameter shows a 
clear correlation with Mbh, which dominates any 
trend with Mi. 

Motivated by these qualitative results, we fit a 
power law of the form 

\ogJ = A + B log (^^) + C(M, + 23) + 

to all SFoo and r data points in each ugriz band 



separately, keeping B fixed to —0.479 and 0.17, 
respectively, in order to avoid the bias discussed 
in Section 15.11 While there is a lot of scatter in 
the variability parameters for fixed M^, Mbh, and 
Ai?F (see Section H?^ . Eq. [7] describes the mean 
trends well. The best-fit coefficients, averaged over 
the five bands, are reported in the first and sixth 
rows of Table [T] The best-fit coefficients when si- 
multaneously fitting all bands are consistent with 
these averages. The reported error bars are com- 
puted from the variation of the best-fit parameters 
over the five bands. 

The dependence of r on redshift and Mi is only 
marginally detected, and can be attributed to the 
L-z degeneracy: as redshift increases, the best-fit 
r decreases, while the coefficient for Mi indicates 
that as luminosity increases, r increases. When 
E is fixed to zero, the dependence on Mi has low 
significance, while the correlation with Mbh re- 
mains the same (seventh row of Tabled]). For il- 
lustration, when E is fixed to 1 [so t oc (1 -I- z), 
third and eighth rows], a spurious dependence on 
Mi emerges. Therefore, the controlling variable 
for the characteristic time scale is clearly Mbh, 
suggesting that more massive black holes vary on 
longer time scales. Similarly, with E fixed to 
for SFoo (second row), the C {Mi) and D {Mbh) 
coefficients remain largely unchanged, confirming 
that there is no significant correlation between am- 
plitude and redshift. Therefore, we force the fi- 
nal adopted model to have no redshift dependence 
{E = 0). Table [1] also provides fits for when D 
is fixed to zero so that SFoo and t can be esti- 
mated in the absence of black hole mass informa- 
tion (fourth and ninth rows). 

Of the analyzed parameters, the probable er- 
rors for black hole mass estimates are the largest: 
of order 0.2 - 0.4 dex (Marconi et al. 2008; Vester- 
gaard & Peterson 2006). When ignored, these 
large statistical uncertainties result in underesti- 
mated values for D (Kelly 2007). With an as- 
sumed random uncertainty of 0.2 dex in the black 
hole mass measurements, we find the best-fit co- 
efficients reported in the fifth and tenth rows of 
Table [1] (continuing to assume that E = Q). It 
can be seen that including mass uncertainties of 
0.2 dex increases D from ^0.1 to ^0.2 for both t 
and SFoo (a bias of about 2-3 formal statistical fit- 
ting errors). The coefficients in the fifth and tenth 
rows represent our best-fit model for the variabil- 
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ity parameters, and we assume these values in the 
remaining Sections. 

The best-fit model shows a smaUcr slope for 
the correlation between r and Mbh of 0.21 ± 0.07 
than the value of 1.0 ± 0.4 found by KBS09. The 
two results are still marginally statistically consis- 
tent given the large uncertainties in KBS09. More- 
over, their sample contains relatively lower mass 
and luminosity quasars than those analyzed in this 
study, and this difference might result in addi- 
tional biases. KozlO also noted that the KBS09 re- 
sults may be affected by contamination from host 
galaxy emission, a problem which will be far smal- 
ler for the generally more luminous SDSS quasars. 

5.3. The Eddington Ratio as the Driver of 
Variability? 

Since the Eddington ratio, L/LEdd, is depen- 
dent on luminosity and black hole mass, the trends 
for SFoo in Table [T] might be explained if SFoo is 
simply driven by L/LEdd- To test this, we es- 
timated L/LEdd as the ratio of the bolometric 
luminosities from Shen et al. (2008) to the Ed- 
dington luminosity, LEdd = 1-5 x lQr'^{MBH / Mq) 
erg/s. Figure [15] shows L/LEdd as a function of 
Mbh and A/^, and demonstrates that lines of con- 
stant L/LEdd are similar in slope to those of con- 
stant SFa^ (see Figure [T4|). The median SF^^ ver- 
sus the median L/LEdd for each bin is shown in 
the right panel, where we find a power law slope 
of —0.23 ± 0.03. This significant anti-correlation 
was also found by Wilhite et al. (2008), Bauer et 
al. (2009), and Ai et al. (2010). KBS09 did not 
report such an anti-correlation; however, they did 
not compare L/LEdd with SF^a but rather with 
SFao/y/r, the driving amplitude of short term 
variability. If L/LEdd is the sole driver of the 
quasar variability amplitude SFoo, we would ex- 
pect that the coefficients for Mi and Mbh are 
related as D = 2.5C. However, we find D = 
(1.37±0.23)C, 4.7cr away from the presumed slope 
of 2.5. This suggests an additional source of vari- 
ability, such as a dependence on luminosity or on 
Mbh in addition to the dependence on L/LEdd- 
Moreover, if we ignore this additional source, the 
result that SFoo depends on L/LEdd supports a 
scenario where the amplitude of the optical vari- 
ability is determined by the accretion rate (see dis- 
cussion in Wilhite et al. 2008). 

If the Eddington ratio is a proxy for the quasar's 



age (e.g.. Martini & Schneider 2003), then a lower 
Eddington ratio, and thus a larger amplitude of 
variability, could indicate a dwindling fuel supply 
and a more variable rate at which it is supplied 
to the black hole. However, it is unlikely that the 
observed variability is due to fluctuations in the 
external fuel supply to the disk, because the fluc- 
tuation time scales of 10 - 10,000 days are very 
short compared to the viscous time scales of 10^ 
- 10^ days that should control large scale changes 
in the accretion rate. Moreover, these fluctuations 
in the fueling rate will also be smoothed out and 
damped as they travel inward, and will likely have 
effectively disappeared by the time they reach the 
optical emitting region (e.g., see discussion in Chu- 
razov et al. 2001). Instead, the origin of the fluc- 
tuations is probably more local. 

Another possibility is that the dependence on 
L/LEdd may simply be a reflection of the depen- 
dence on wavelength, which in turn depends on 
the disk radius. If a higher L/LEdd means a 
hotter disk, then the optical flux originates at a 
larger radius. Assuming that longer wavelengths 
are emitted further out in the disk, the lower 
SFoo at longer wavelengths would lead to the anti- 
correlation between S'-F'oo and L/LEdd- In thin 
disk theory (e.g., see Frank, King, & Raine 2002), 
the characteristic radius for emission at wave- 
length Xrf scales as R\ cx m'^J^{L / LEddY^^ X^^i/^, 
and the thermal time scale is related to the orbital 
time scale as tth oc torb oc R^^^ /^/Mbh- There- 
fore, under the assumption that r is related to the 
thermal time scale, and that variability at wave- 
length XfiF is dominated by the scale R\, r scales 
with XnF, Mbh, and L as: 

r^M'J^{L/LEddy/^X'RF- (8) 
Since LEdd oc Mbh, this means r scales as 

t^L'/'XIf, (9) 

which does not match the observed scaling of r cx 
^--0.075 ;)^ai7^j0.2i^ If we assume that r is related 

to the viscous time scale, 

rocLy^'X^I'^Ml^ (10) 

still in conflict with the measured values. 

The variability time scale simply does not show 
the strong dependence on Xrf, Mbh, and L ex- 
pected from these simple scalings. Therefore, us- 
ing this naive scaling, we are not able to relate the 



12 



Table 1 

Best-Fit Coefficients for Eq. [7] 



/ 


A 


B (Xrf) 


C (M,) 


D (Mbh) 


E {1 + z) 




-0.57 ±0.01 


-0.479 ± 0.005 


0.117 ± 0.009 


0.11 ± 0.02 


0.07 ± 0.05 


SFoo 


-0.56 ±0.01 


-0.479 ± 0.005 


0.111 ± 0.005 


0.11 ± 0.02 


= 


SFac 


-0.760 ± 0.009 


-0.479 ± 0.005 


0.193 ± 0.006 


0.12 ± 0.02 


= 1 


SFoc 


-0.618 ±0.007 


-0.479 ± 0.005 


0.090 ± 0.003=^ 


= 


= 




-0.51 ± 0.02 


-0.479 ±0.005 


0.131 ± 0.008'^ 


0.18 ± 0.03 


= 


T 


2.4 ± 0.2 


0.17 ± 0.02 


-0.05 ± 0.03 


0.12 ± 0.04 


-0.7± 0.5 


T 


2.3 ±0.1 


0.17 ± 0.02 


0.01 ± 0.03 


0.12 ± 0.04 


= 


T 


2.1 ±0.1 


0.17 ±0.02 


0.09 ± 0.03 


0.13 ± 0.04 


= 1 


T 


2.2 ±0.1 


0.17 ± 0.02 


-0.01 ± 0.02 


= 


= 



2.4 ±0.2 0.17 ±0.02 0.03 ± 0.04 0.21 ± 0.07 =0 



^Bascd on K-corrected magnitudes. Without the K-correction, this cocfhcent changes to 
0.079 ± 0.003. 

''Based on K-corrected magnitudes. Without the K-correction, this cocfhcent changes to 
0.113 ± 0.006. 

"^Measurement errors in Mbh of 0.2 dex have been included in the fitting. These coefli- 
cients are recommended when applying the model. 

Note. — In each row, the B coeflicient was determined and fixed before fitting a multiple 
regression in all other parameters (see Section lS.lll . The cosmology used for determining Mi 
is Qm = 0.30, = 0.70, and h = 0.70 (Schneider et al. 2007), whereas that used in the 
Mbh estimates is Qm = 0.26, Ha = 0.74, and h = 0.71 (Shen et al. 2008). This difference 
should only have a 1% effect on the best-fit coefficients. 
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observed r to either a thermal or viscous time scale 
of the radius associated with the wavelength of the 
variabihty. However, in reality there is a range of 
radii, which corresponds to a range in time scales, 
contributing to the observed flux in each band, and 
this will cause some degree of smoothing. Also, the 
radial regions might overlap for each band, caus- 
ing a single radius to contribute flux in multiple 
bandpasses, and this might attenuate some of the 
dependencies on Xrf, Mbh, and L. Finally, our 
observed power law indices may be biased due to 
the uncertainty in the bolometric correction. 

5.4. A "Recipe" for Generating Mock 
Light Curves 

Given the success of the damped random walk 
model in explaining a large body of observations, 
it represents a simple quantitative framework for 
generating mock quasar light curves. The optical 
light curves can be simulated for quasars at dif- 
ferent redshifts and for a wide range of luminosity 
and black hole mass, which provides the basis for 
quantitative modeling and optimization of quasar 
variability surveys. The model presented has al- 
ready been implemented as a part of the simula- 
tion effort in support of the Large Synoptic Survey 
Telescope (Ivezic et al. 2008). 

At a given redshift, a simulated value for ab- 
solute magnitude can be readily drawn from an 
adopted luminosity function (e.g., Croom et al. 
2009 and references therein). An estimate for the 
black hole mass is also required in order to ap- 
ply our model, and it can be generated using the 
adopted absolute magnitude as follows. Accord- 
ing to the Shen et al. (2008) results, the quasar 
luminosity and black hole mass are strongly cor- 
related (see the bottom-right panel of Figure [T^ . 
Using Shen et al. (2008) values, we quantify this 
correlation in Figure [T51 At a given M^, the Mbh 
distribution has a finite width due to the distribu- 
tion of Eddington ratios and measurement errors. 
Both effects can be well-described using a Gaus- 
sian distribution: 

p{\0gMBH\M,) = 
1 r (l0gMBff-/i)'l 



where = 2.0 - 0.27M, and cr = 0.58 + O.OllAf, 
(the black hole mass is expressed in solar mass 
units). 



With adopted values for Mi and Mbh-, the vari- 
ability amplitude and the characteristic time scale 
can be estimated using Eq. [T] We note that choos- 
ing a time scale from an adopted value for SFao 
by utilizing the correlation seen in Figure [3] (and 
drawing from a Gaussian distribution similar to 
Eq. [TT]) is not as accurate as using Eq. [7] and the 
adopted values of Mbh and \b,f- Finally, given 
a mean magnitude and wavelength, a quasar light 
curve can be easily generated using Eq. [5j 

6. Variability Properties of Radio- and X- 
ray detected Quasars 

Previous studies, based mostly on two-epoch 
datasets, found that radio-loud quasars are margi- 
nally more variable than radio-quiet quasars for 
rest-frame time lags in the range 50-400 days (e.g., 
VB04 and references therein). It was also found 
that X-ray-detected quasars are significantly more 
variable than X-ray-undetected quasars at rest- 
frame time lags up to 250 days. 

The large size and the availability of light curves 
for our sample allows us to revisit these results 
with more statistical power. In addition, the 
damped random walk model and our best-fit re- 
sults from Table [T] provide a convenient method 
to account for various selection effects. For ex- 
ample, it is possible that subsamples selected by 
various means, such as the requirement of radio 
or X-ray detections, have different distributions 
of luminosity and black hole mass. In this case, 
they would display different variability behavior 
not because of intrinsic differences in the variabil- 
ity mechanism, but rather because of the trends 
captured by Eq. [71 A comparison of the ratio of 
observed and predicted variability parameters for 
any subsample and the full sample will automati- 
cally take into account these sampling effects, and 
this is the main statistical method we use in this 
section. We first discuss a radio subsample, and 
then analyze a sample of quasars with X-ray de- 
tections. 

6.1. Radio-detected Subsample 

We use the unified radio catalog of Kimball & 
Ivezic (2008) to access the FIRST and NVSS data 
(20 cm continuum data) for our objects. We refer 
the reader to this paper and references therein for 
details about these radio surveys and object asso- 
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ciation. Figure [17] shows the radio-detected frac- 
tion of quasars in S82 as a function of magnitude. 
The overall fraction (5%) is considerably lower 
than that in the Northern galactic cap footprint 
of the SDSS for the magnitude range 19 < « < 20 
(White et al. 1997; Ivezic et al. 2002). This is due 
to a difference in targeting algorithms between the 
two surveys, and results in nearly 300 quasars with 
both optical light curves and radio data. When 
available, we use NVSS values for radio flux, but 
otherwise we adopt the FIRST integrated flux. 
The radio-loudness parameter Rg is calculated ex- 
actly as the Ri in Ivezic et al. (2002), but using g 
magnitudes instead. We assume a spectral index 
of —0.5 for the optical and for the radio (because 
this sample is dominated by the core radio emis- 
sion; see Kimball & Ivezic 2008) when computing 
the K-correction. 

Table [2] compares the mean properties for vari- 
ous subsamples detected in the radio to the radio- 
undetected subsample (essentially the full sample 
due to small fraction of radio-detected objects). 
For each subsample, the mean t and SF^o are 
listed, as well as the mean ratio of the observed 
values to those predicted using Eq.[7]and the mea- 
sured values of \rf, Mi, and Mbh- The mean 
values are computed iteratively with ±3(t outliers 
excluded. Errors are reported as the (clipped) rms 
divided by \/N — 1, where N is the number of data 
points in the distribution. Numbers of objects are 
lower for columns involving model quantities be- 
cause black hole mass estimates are not available 
for all objects. 

Table [2] shows that only the variability ampli- 
tude of the radio-loudest quasar subsample (Rg > 
3) is significantly different (> 3a deviation) from 
the behavior of the full sample. The radio-loudest 
quasars have systematically larger variability am- 
plitudes, when corrected for trends described by 
Eq. [71 by about 30%, compared to the full sample 
dominated by radio-quiet objects, in agreement 
with VB04. 

6.2. X-ray Detected Subsample 

For the analysis of variability properties of 
quasars detected at X-ray wavelengths, we use 
data from the ROSAT All-Sky Survey (RASS; Vo- 
ges et al. 1999). The X-ray subsample consists 
of 82 quasars with RASS full-band count rates 
greater than 10^'^ ct/s, taken from the Schnei- 



der et al. (2007) catalog. As can be seen from 
Table [H the variability properties of this subsam- 
ple are statistically indistinguishable from the full 
sample. VB04 detected a significant increase in 
structure function at rest time lags below 250 days 
for their X-ray subsample, but on long time scales 
there was no significant difference. It is plausible 
that their result for short time scales was infiu- 
enced by small sample size and t effects discussed 
in Section [331 

7. Summary and Conclusions 

We have used the damped random walk model 
of KBS09 and KozlO to model the optical/UV 
variability of ^9000 SDSS Stripe 82 quasars with 
the ugriz light curves. The dataset includes 2.7 
million photometric measurements collected over 
10 years. We confirm that this is a good model 
of quasar variability, and quantify the dependence 
of two variability parameters, the long-term rms 
variability SFoo and the damping time scale t, on 
physical parameters such as wavelength, luminos- 
ity, black hole mass and Eddington ratio. Our 
main results are the following: 

1. A stochastic process with an exponential co- 
variance function characterized by an am- 
plitude and time scale provides a good fit 
to observed quasar light curves, as shown 
by KBS09 and KozlO, using smaller sam- 
ples with less wavelength coverage, but bet- 
ter time sampling. 

2. The long-term rms variability, SFoo, has a 
mode at ~0.2 mag and characteristic time 
scales, T, are roughly 200 days in the rest 
frame, as found previously by KBS09 and 
KozlO. These time scales are consistent with 
thermal time scales, but simple accretion 
disk models fail to reproduce the observed 
scaling of r with physical parameters. 

3. Quasars with similar physical parameters 
can have different characteristic time scales 
for variability. It is now clear that the dis- 
tribution of T accounts for most of the scat- 
ter in the structure function on short time 
scales for quasars with similar luminosity, 
rest wavelength, and time lag, which ex- 
plains the puzzling results from MacLeod 
et al. (2008). Results from fitting a power 
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Table 2 

Mean Variability Properties of Radio and X-ray Subsamples 





< log T> 


N 


<log(T/Tp)> 


N 


<logSFoo> 


N 


<\oe{SF^/SF^_p)> 


N 


no radio 


2.305 ± 0.005 


6989 


-0.003 ± 0.007 


4467 


-0.634 ± 0.003 


7067 


0.018 ± 0.003 


4508 


no radio i < 19 


2.327 ± 0.013 


1519 


-0.013 ± 0.013 


1279 


-0.721 ± 0.006 


1527 


0.003 ± 0.006 


1287 


radio 


2.26 ± 0.03 


277 


-0.07 ± 0.05 


154 


-0.70 ± 0.02 


283 


-0.03 ± 0.02 


154 


radio i < 19 


2.24 ± 0.06 


105 


-0.08 ± 0.06 


73 


-0.81 ± 0.03 


108 


-0.08 ± 0.03 


74 


Rg > 3 


2.30 ± 0.05 


135 


-0.01 ± 0.07 


68 


-0.57 ± 0.02 


137 


0.11 ± 0.03 


68 


-R, < 3 


2.28 ± 0.03 


211 


0.00 ± 0.04 


133 


-0.73 ± 0.02 


219 


-0.06 ± 0.02 


140 


fig < 2 


1.99 ± 0.14 


30 


-0.21 ± 0.14 


19 


-0.88 ± 0.07 


30 


-0.12 ± 0.06 


19 


resolved'^ 


2.29 ± 0.05 


111 


0.00 ± 0.06 


63 


-0.63 ± 0.02 


111 


0.04 ± 0.03 


63 


unresolved^ 


2.24 ± 0.04 


167 


-0.13 ± 0.07 


91 


-0.73 ± 0.03 


172 


-0.06 ± 0.03 


91 


x-ray 


2.41 ± 0.05 


81 


0.03 ± 0.06 


59 


-0.60 ± 0.03 


82 


0.01 ± 0.03 


58 


no x-ray 


2.307 ± 0.005 


6950 


-0.004 ± 0.007 


4559 


-0.638 ± 0.003 


7020 


0.017 ± 0.003 


4598 





^The morphological radio classes are defined using the integrated and peak FIRST fluxes as in Kimball & Ivezic (2008). 



Note. — r and SFao are the observed time scales and asymptotic amplitudes of optical variability (Section |4)l, while Tp 
and SFcx>,p refer to those predicted from Eq. O using the coefficients in the fifth and the tenth rov^s of Table [T] Errors 
are reported as the rms divided by ^^N — 1, where N is the number of data points, listed after each column. Numbers of 
objects are lower for columns involving model estimates because black hole mass estimates are not available for all objects. 
Only the variability amplitude for the radio-loudest quasar subsamplc {Rg > 3) is significantly different (> 3cr deviation) 
from the behavior of the full sample. 



law to observed ensemble structure functions 
should be interpreted with caution. 

4. The variability time scale is correlated with 
the long-term rms variability with a slope 
of 1.30 ± 0.01 dex/dex. Quasars that have 
large long-term variability amplitudes gener- 
ally vary on longer characteristic time scales. 
The amplitude of short term variations is 
also correlated with r. This conclusion is 
unaffected by any time sampling issues in 
the S82 dataset. 

5. The damped random walk corresponds to 
a PSD proportional to 1//^ at frequencies 
/ > (27rT)^^, flattening to a constant at 
lower frequencies. At large /, the data are 
in great agreement with PSD cx 1//^. In 
terms of the structure function, this means 
that SF (X (At)^/^. Whereas previous anal- 
yses of the SF obtained a power law slope 
of (3 = 0.3, here we demonstrated that this 
would be a consequence of fitting the data 
around the "knee" (turn-over) of the SF. 
Our constraints for small / are much weaker. 
As discussed in Section 14.41 due to a lack 
of sufficient long-time scale information, we 
are unable to distinguish between a 1//° or 
a 1// PSD at frequencies / < (27rT)^^ us- 



ing the data and computational technique 
described here. 

6. The rest-frame variability parameters show 
a negligible trend with redshift, suggesting 
that they are intrinsic to the quasars, and 
these properties do not evolve over cosmic 
time for fixed physical parameters of the 
quasar {AIbh, Mi, and Xrf)- 

7. For fixed luminosity and black hole mass, 
T increases with increasing rest-frame wave- 
length with a power law index of 0.17, and 
SFoo decreases with a power law index of 
—0.48. The latter result is similar to pre- 
vious findings (e.g., MacLeod et al. 2008, 
VB04). KozlO also observed that the vari- 
ability increases to shorter A, but they kept 
r fixed in their fits. If wavelength is a proxy 
for radius in the accretion disk, this implies 
that the characteristic time scales are longer 
and the variability amplitudes are smaller in 
the outer regions than in the inner regions. 

8. The long-term variability SF^o is strongly 
anti-correlated with luminosity as found in 
previous studies such as VB04, Wilhite et 
al. (2008), and references therein. By study- 
ing the median SF,^ in the plane of ab- 
solute magnitude and black hole mass, we 
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can separate the anti-correlation of ampli- 
tude with luminosity from the positive cor- 
relation with black hole mass. As suggested 
in Wilhitc ct al. (2008), these trends may 
be largely explained if the amplitude of vari- 
ability is tied to changes in the accretion rate 
in the disk, and is simply related to the Ed- 
dington ratio. However, despite the strong 
anti-correlation between SFoc and L/Lsdd 
(which accounts for most of the dependence 
on Mi and Mbh), the exact dependence 
with Mi and Mbh is not consistent with 
L/LEdd as the sole driver of quasar variabil- 
ity. 

9. The damping time scale r appears to be 
nearly independent of luminosity and corre- 
lated with Mbh with a power law index of 
0.21 ± 0.07. The mild discrepancy with the 
KBS09 result (1.0±0.4) may be due to the 
different range of sampled luminosity and 
black hole mass, as well as contamination 
by host galaxy emission in many of the very 
low luminosity systems they consider (see 
KozlO). 

10. While the mean variability parameters can 
be related to physical parameters, for fixed 
values of Mj, Xrf, and Mbh, there is still 
a large scatter around the mean values, sim- 
ilar to the variance of the observed distri- 
butions. Some of that scatter can be at- 
tributed to measurement and fitting errors 
(~60% for T and ~70% for SF^), but there 
is enough evidence for residual stochastic na- 
ture of quasar variability. Therefore, it can- 
not be assumed that quasars with similar 
Mil Xrf, and Mbh will necessarily have 
similar variability properties. 

11. The radio-loudest quasars have systemati- 
cally larger variability amplitude by about 
30%, while the distribution of their charac- 
teristic time scale is indistinguishable from 
that of the full sample. There are no statis- 
tically robust differences in the characteris- 
tic time scale and variability amplitude be- 
tween the full sample and a small subsample 
of quasars detected by ROSAT. 

With this paper, and results from KBS09 and 
KozlO, the ability of the damped random walk 



model to quantitatively describe quasar variabil- 
ity is well-established. As emphasized by KozlO, 
this means that variability studies can become 
fully quantitative because the entire process of 
identifying and assigning parameters to quasars 
can be simulated to allow estimates of complete- 
ness and parameter biases. In particular, an im- 
portant next step is to determine the variability 
equivalents of luminosity functions, i.e., the in- 
trinsic distributions of the variability parameters. 
While our results represent a good first step in 
this direction, we caution that a non-negligible 
fraction of the S82 quasars have indeterminately 
long time scales. If, following KBS09, we iden- 
tify the characteristic time scale with the thermal 
time scale, then the next question is whether there 
are additional time scales (such as the dynami- 
cal or viscous time scale), or sources of variabil- 
ity. For example, Blackburne & Kochanek (2010) 
recently found evidence for changes in disk size 
with changes in luminosity using gravitational mi- 
crolensing. Such searches for additional sources of 
variability will likely require better sampled light 
curves and over a longer time scale. 

The prospect of advancing these studies of 
quasar variability now faces a bottleneck. The 
S82 quasars have the advantage of sample size, 
wavelength coverage, and spectroscopy, but the 
light curves have poor sampling and modest over- 
all lengths, leading to significant problems for 
accurately estimating r when it is long. The 
quasars behind the Magellanic Clouds (Kozlowski 
& Kochanek 2009) are a smaller sample without 
spectroscopic confirmation, but have superb, long 
term light curves that continue to be extended be- 
cause of the continuing microlensing projects. Im- 
proving on our present results in the short term de- 
pends on either reviving the monitoring of S82 or 
spectroscopically confirming the Magellanic Cloud 
quasars. 

Resuming the monitoring of Stripe 82 is chal- 
lenging with the decommissioning of the SDSS 
imaging system. The best short-term prospects 
are the Pan-STARRS project (Kaiser et al. 2002), 
the Palomar-QUEST project (Schweitzer et al. 2006), 
or using the DECam being built for the Dark 
Energy Survey (Honscheid et al. 2008). Since 
the challenge is to constrain long time scales, the 
presence of a multiple- year gap is mainly a com- 
plication for ensuring that problems in matching 
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photometric bands are not interpreted as a form 
of long-term variability. Obtaining spectra of the 
Magellanic Cloud quasars is in some ways easier, 
because the quasar magnitudes and densities are 
well-suited to the AAOmega fiber spectrograph 
on the A AT and, to a lesser degree, the IMACS 
spectrograph on Magellan. 

The more general. Northern monitoring projects 
such as Pan-STARRS and Palomar-QUEST will 
slowly build light curves for essentially all the 
SDSS quasars, but at present their cadences are 
not ideal (see Schmidt et al. 2010) and it will take 
nearly a decade to build the long duration light 
curves needed for the analysis. In the long term, 
observations will be significantly improved with 
the advent of next-generation sky surveys. Most 
notably, the Large Synoptic Survey Telescope 
(LSST, Ivezic et al. 2008) will obtain accurate, 
well-sampled light curves for millions of AGN. 
The observed distribution of rest-frame character- 
istic time scales for S82 quasars spans the range 
from about 10 days to 1000 days (c.f. Figure [3]). 
To probe the time scales as short as O.lr, and 
assuming a characteristic redshift of 2, the light 
curves should be sampled every 3 days in the ob- 
server's frame, which is in good agreement with 
the baseline cadence of LSST. With a 10 year-long 
survey, the length of the light curves will be in the 
range (1-200)t. A combination of the SDSS, Pan- 
STARRS, DBS, and LSST data for -10,000 Stripe 
82 quasars would span well over two decades, with 
multi-band photometry obtained for hundreds of 
epochs, and would represent the best sample to 
date for studying the optical continuum variabil- 
ity of quasars. In particular, such a dataset would 
enable a robust measurement of the low-frequency 
behavior of their PSD (c.f. point 5 above). For il- 
lustration, the LSST photometric errors in the r 
band will be < 0.02 mag for r < 22, and there 
are roughly 2-3 million AGN with r < 22 in the 
20,000 sq. deg. covered by the main LSST survey 
(see Table 10.2 in the LSST Science Book; Abell et 
al. 2009). Each of these objects will be observed 
about 1000 times, yielding a database of over 2 
billion photometric measurements. This data set, 
roughly a thousand times larger than that ana- 
lyzed here, will enable a significant improvement 
in our understanding of quasar variability. 
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Appendix: Search for Periodic Light Curves 

Although a stochastic process has proven to be an accurate statistical description of quasar light curves, 
any discovery of periodic behavior (even of a single source) would have interesting physical implications. 
However, a periodogram can only be used as an indicator of significant periodicity in a signal as compared 
to pure white noise (i.e., having no signal at all). Since quasars are genuinely variable, as described by a 
damped random walk (i.e., a red noise process; see KBS09), one can only evaluate the significance of the 
periodicity when also allowing for the signal covariance that is also present (see Markwardt et al. 2009; Gotz 
et al. 2009; Cenko et al. 2009, and references therein). Nevertheless, as a quick test for outstanding cases of 
periodicity among the S82 light curves, we analyzed 8,863 light curves for evidence of periodicity using the 
Lomb-Scargle periodogram (Lomb 1976; Scargle 1982). The threshold for considering the strongest peak in 
the periodogram as candidate evidence for periodicity was set following Horne & Baliunas (1986), with an 
adopted false alarm probability of 0.05. However, this is the threshold for ruling out white noise in favor 
of periodicity, and therefore provides no information on how a periodic description compares to a colored 
noise process such as a damped random walk. To determine the latter, a different threshold is needed (see 
KozlO). When adopting the threshold for ruling out white noise in favor of periodic variability, we identify 
88 light curves as good candidates. 

A close inspection of the period distribution for the full sample and the selected 88 candidates shows 
important differences: while the full sample displays a fairly flat distribution ranging from 100 days to values 
exceeding 10,000 days, the period distribution for 88 candidates is bimodal. The first peak with 22 objects 
corresponds to aliasing at roughly one-year sampling cadence, while the second peak is centered on periods 
of about 6-7 years, similar to the total length of observations. We have visually inspected the light curves 
and phased light curves for all 88 candidates. It turns out that candidates with proposed periods of the 
order one year have light curves consistent with aliasing, while those with longer periods typically have only 
one observed "oscillation" that might not be used as robust evidence for periodicity (six examples are shown 
in Figure [TS]) . Therefore, our search for periodic light curves in the S82 quasar sample has not yielded 
any convincing cases. Again, it is important to note that in cases such as in Figure I18[ the periodogram 
indicates that these sources with long r are likely to be periodic, not because they show true periodicity, but 
because it computes the likelihood relative to the wrong null hypothesis (white noise rather than colored 
noise). Moreover, the fact that <5% of the light curves exceed the threshold for periodicity being a better 
description than pure white noise is a further indication that the periodogram is not a very powerful statistic 
for poorly sampled light curves, since it is clear from the previous Sections that quasar light curves are 
not white noise, but rather are well described by a damped random walk. In fact, the single "oscillation" 
observed for the best candidate long period objects is entirely expected from colored noise processes such as 
a damped random walk (see right panels of Figure [T5)) . Therefore, this analysis is further evidence that a 
damped random walk is a good description of quasar light curves. 
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Fig. 1. — The top panel shows a segment of a simulated light curve with t — 20 days and SFoc — 0.10 
mag. In the bottom panel, triangles, squares, and crosses represent the SF computed for the first, second, 
and third sections of the total light curve, respectively (each has a length of 500t). The increased scatter in 
SF{At) at large time lags is due to the finite length of each light curve section - the scatter decreases when 
using the entire light curve length of ISOOr (as shown by the solid line). The short-dashed line is a power 
law fit, SF{At) oc At^ with P = 0.3, to the data points with 0.15 < At/r < 3. The long-dashed line is the 
true function SF = 5Foo[l - e-"^*/^]^^^. 
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Fig. 2. — Initial light curve selection. Top-left: the distribution of the number of observations per light 
curve for the r (solid) and u (dashed) bands. Top-right: distribution oi ALnoise- We define light curves with 
ALnoise < 2 to be more consistent with uncorrelated noise rather than our model. Bottom-left: distribution 
of ALoo', light curves with ALoo < 0.05 likely have run-away time scales. In this panel and the previous 
panel, the x-axes are truncated at 0.5 and 5 respectively, but the histograms continue to greater values. 
Bottom-right: distribution of psr degrees of freedom (Ndof) for the damped random walk model (solid 
line) . The expected Gaussian distribution based on Ndof is also shown (dashed) . The hashed region in each 
panel shows the values rejected from our final sample. 
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Fig. 3. — Top: Distributions of the best-fit rest-frame time scale r (left) and long-term structure function 
•S'-F'oo (right) for the S82 quasars (33, 112 light curves). The filled histograms show the fits for quasars with 
absolute i-band magnitude in the range —27 < Mi < —26 (the median Mi for the whole sample is —25.46). 
Bottom: Relationship between r and SF^o- A power law shown by the dashed line is fit to the medians 
(open white circles), with the slope listed in the top-right corner. Contours show regions containing 90%, 
70%, 50%, and 20% of the total number of points. The star symbols label the four regions in parameter 
space from which the sample light curves in Fig. 2] (labeled A through D) were chosen. The observed-frame 
r distribution lacks many of the short time scales observed in KozlO. This is likely due to either the better 
time sampling of the OGLE light curves or stellar contamination in their sample. 
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Fig. 4. — Each panel shows a sample light curve from a certain region of r, SFoo space (the four regions A-D 
are indicated with stars in the bottom panel of Figure [S]) . Data points with error bars show the observed 
S82 data. Solid lines show the weighted average of the damped random walk model light curves that are 
consistent with the data (see Section[2]). Dotted lines show the ilcr range of these possible stochastic models 
about the average. 
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Fig. 5. — Test of biases in best-fit parameters. Left: The input variability parameters r and SFoa were drawn 
from a uniform distribution in log r and log SFoo limited by the dotted rectangle. The contours show the 
output distribution after applying the cuts described in Section [3.11 and the dashed line is a linear regression 
between the output parameters (with the slope listed on the panel). Solid, blue circles show the objects 
which do not satisfy ALqo > 0.05 (these make up 21% of the starting sample, where 13% are saturated at 
T = lO'^ days). The small black dots show the input values for these rejected points. Open, red circles show 
the 3% that do not satisfy ALnoise > 2. Right: Relationship between Tout /Tin and S Font /S Fin- The slope 
of the cleaned output distribution is listed. The sharp edge is due to the saturation limit of r = 10^ days. 
The contours show the 90%, 70%, 50%, and 20% levels. 
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Fig. 6. — Top panel: Data points with error bars show the observed hght curve for a quasar in our final 
sample, the solid line shows the weighted average of all consistent model light curves, and the dotted lines 
show the ±1(7 range of possible stochastic models (see Section [2]). Bottom panel: A "regenerated" light 
curve using Eq. [5] and the estimated values of t and SFoo listed in the top panel. Due to the poor time 
sampling of the light curve, the best fit to the "regenerated" light curve has a run-away time scale (listed at 
the bottom), with AL^ = 0.003. 
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Fig. 7. — Results from "regenerating" light curves derived from the observed t and SF^o (see Section l4?2l) . 
The solid, black lines in the top panels show the distribution of input SFoo (left) and r (right) normalized 
by the median input value, and the ratio of output to input values are shown by the red, dashed lines. The 
bottom panels show the scatter in the output-to-input ratios (left), and the relationship between the output 
parameters (right). Listed in each bottom panel is the slope of a linear fit (dashed line) to the median values 
(open circles). Contours show the 90%, 70%, 50%, and 20% levels. 
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Fig. 8. — The symbols with error bars show the distribution of S Fobs / S Fmodei froni Figure 2 in MacLeod et 
aL (2008), where SFobs is the observed structure function for individual S82 light curves for a fixed observed- 
frame time lag of 1 year, and SFmodei = 1(1 + 0.03Afi)(AiflF/AflF)° '*^- The distribution from the ensemble 
analysis by 104 (based on two-epoch SDSS data) is shown by the dashed curve, scaled by a factor of 0.2. The 
dotted histogram is the distribution of ■\/mcdian(T)/T from this work, where r is the characteristic time scale 
for each light curve. Its agreement with the S Fobs I S Fmodei distribution for individual light curves shows 
that most of the scatter pointed out by MacLeod et al. (2008) is due to a finite width of the r distribution. 
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Fig. 9. — First three panels: Simulated light curves using r = 5.8 days, rms = 0.124 mag, and a PSD 
described by PSD oc with a = —2 for / > (27rr)~^, and a = —1 {red; top-left), a ~ {blue; top-right), 
or a = —1.9 {green; bottom-left) for / < {2ttt)^^ . The lines in the last panel show each PSD from which the 
light curves are generated, with the black dotted line indicating the break frequency. As a check, the PSD 
was then computed for the simulated light curves, shown with colored dots. The y-axis in the bottom-right 
panel is in arbitrary units. 
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Fig. 10. — In the top panels, the filled histogram shows the ("input") distributions of t and = /^dof ~ X^/^dof for 

r-band S82 light curves (~7000 total). This distribution of r, along with the r-band rms values, are used to generate ^TOOO 
realizations of noise processes with a = —2 for / > (27rr)~^ and a = 0, a = —1, or a = —1.9 for / < {2ttt)~^ . The dashed 
blue, solid red, and dotted lines show the best-fit ("output") distributions when modeling these three processes, respectively, 
as a damped random walk (DRW). Here, the simulated observations are spaced every 5 days over 10 years with typical errors 
of 0.01 mag. The output a, SFoo, and t are compared to the input values in the next three panels. In the last panel, the 
correlation between the output log{SFoo) and logij) is shown using a linear fit to the median values {open circles), with the 
slopes listed in the legend for each case. Contours show regions containing 90% of the data points. 
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Fig. 11. — As in Figure [TOl except with the Stripe 82 time samphng and photometric accuracy imposed on 
the simulated light curves. 
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Fig. 12. — The top panels show the distribution of quasars in redshift (left) and i magnitude (right) on a 
grid of Mbh vs. Mi. The bottom panels show the number of data points in bins of Mbh (or redshift) vs. 
Mi. 
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Fig. 13. — Distribution of the rest-frame time scale r {top panel) and long-term structure function SFoa 
{bottom panel) as a function of rest-frame wavelength Xrf- The different contours show the 70% and 30% 
levels for the u, g, r, i, and z (from left to right) bands. The best-fit power law for each subsample is shown 
with the same line style as the contours. The thick red line connecting the left and right axes shows the 
median power law slope derived by fitting individual quasars (see text), and has a value of i? = 0.17 and 
—0.479 for T and SF^, respectively. The slope within each band from fitting ensembles of quasars can be 
quite different from the overall slope because of the L-z degeneracy (see text). 
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Fig. 14. — Top panels: The long-term rms variability SFoo (left) and characteristic time scale r (right) are 
shown as colors on a grid of redshift and absolute i-band magnitude Mi. The structure function parameters 
are normalized to a fixed rest wavelength using the fitted power law dependencies of (A_r_f/4000A)^ with 
B = —0.479 and 0.17 for SFoo and r, respectively. The lines of constant variability (dashed) show that 
S'-F'oo is independent of redshift. Bottom panels: As in the top panels but with black hole mass Mbh on the 
X-axis. 
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Fig. 15. — Left: The Eddington ratio for S82 quasars (estimated using masses and bolomctric luminosities 
from Shen et al. 2008) is shown as colors on a grid of Mi vs. Mbh, with dashed lines of constant L/LEdd 
over-plotted. Right: Long-term rms variability (corrected for wavelength dependence) is shown as a function 
of L/LEdd (open circles are medians in each bin). The slope of the linear fit to the medians is listed on the 
panel. 
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Fig. 16. — Determination of p{MBH\Mi) (see Eg. [TT|). The histograms in the top three panels show the 
distribution of black hole masses for S82 quasars with absolute magnitudes within 0.5 mags of the value listed 
in the upper right corners. In the bottom panel, the mean value of log (Mbh/Mq) from each histogram is 
plotted against Mi. 
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Fig. 17. — Top: The number of FIRST radio detections (thick Une, multiplied by 10) as compared to the 
total number of quasars (thin line) as a function of i magnitude in the Northern galactic cap (NGC) footprint 
of the SDSS. Middle: As in the top panel but for Stripe 82. Bottom: Fraction of quasars in the NGC (thin 
line) and in S82 (thick line) samples that have radio detections. Between i = 19 and 20, the fraction of 
radio-detected quasars is considerably lower in S82 than in the NGC. 
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Fig. 18. — Six examples of candidate light curves showing significant periodicity as compared to pure noise (see Appendix). 
The phased light curves are shown on the left, where the phase is the fractional part of the ratio of time to the period (i.e., 
the phase zero-point is arbitrary). The best-fit period is listed on the bottom of each panel. The full light curves are shown on 
the right, where the solid lines show the weighted average of all consistent damped random walk model light curves, and the 
dotted lines show the ztlcr range of possible stochastic models (see Section The best-fit damped random walk parameters 
are listed at the bottom of the right panels, along with the reduced ■ 
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